“Baby Mental Life: Study 2” was conducted on MTurk on 2018-08-04.
Our planned sample was 300 participants, and we anticipated that roughly 80% of recruited participants would pass all of our attention checks, so we initially recruited 378 participants (on the idea that ~80% of 378 ~ 300 participants; note that for administrative purposes we need to recuit participants in batches that were divisible by 9). After filtering out participants who failed at least one of our attention checks, we ended up retaining fewer than 300 participants, so we recruited an additional 16 participants for a total of 394 people recruited. At each stage, we recruited women and men through separate studies, in hopes of acquiring a roughly equal split between genders.
In the end, we ended up with a sample of 304 participants who passed our attention checks, 237 of whom came from unique GPS coordinates.
For this first pass, these data exclude participants where there is another participant with an identical set of GPS coordinates as recorded by Qualtrics.
Each participant assessed children’s mental capacities at 13 target ages between the ages of 0 and 5 years. For each target, they rated 20 mental capacities on a scale from 0 (not at all capable) to 100 (completely capable).
For more details about the study, see our preregistration here.
Here we run some exploratory analyses on these data.
# load required libraries
library(tidyverse)
library(langcog) # source: https://github.com/langcog/langcog-package
library(psych)
library(lme4)
# set theme for ggplots
theme_set(theme_bw())
# run source code (extra home-made functions)
source("./scripts/max_factors_efa.R")
source("./scripts/plot_fun.R")
source("./scripts/reten_fun.R")
source("./scripts/table_fun.R")
source("./scripts/data_prep.R")
NAs introduced by coercionattributes are not identical across measure variables;
they will be droppedJoining, by = "question_qualtrics"
First, I’ll do the factor analysis, and get a list of the top 5 items by factor:
# load in S1 efa in case we need it
efa_S1 <- readRDS("../study 1/s1_efa.rds")
demo_S1 <- read.csv("../study 1/s1_demo.csv")
# conduct S2 efa
efa_S2 <- fa(d_all, nfactors = 4, rotate = "oblimin", fm = "minres",
scores = "tenBerge", impute = "median")
Loading required namespace: GPArotation
table_fun(efa_S2, num_items = 5, pos_abs = "abs")
| capacity | loading |
|---|---|
| Factor 1 | |
| feeling_distressed | 0.82 |
| feeling_overwhelmed | 0.82 |
| feeling_frustrated | 0.78 |
| feeling_helpless | 0.76 |
| feeling_lonely | 0.60 |
| Factor 2 | |
| having_self_control | 0.96 |
| controlling_their_emotions | 0.93 |
| telling_right_from_wrong | 0.93 |
| planning | 0.89 |
| reasoning_about_things | 0.89 |
| Factor 3 | |
| getting_hungry | 0.83 |
| feeling_pain | 0.79 |
| feeling_tired | 0.67 |
| hearing_sounds | 0.58 |
| feeling_physically_uncomfortable | 0.49 |
| Factor 4 | |
| feeling_happy | 0.82 |
| finding_something_funny | 0.81 |
| feeling_excited | 0.79 |
| loving_somebody | 0.64 |
| learning_from_other_people | 0.51 |
Now I’ll use this to define categories of mental capacities:
# get items by factor
factors_S2 <- efa_S2$loadings[] %>%
data.frame() %>%
rownames_to_column("capacity") %>%
gather(factor, loading, -capacity) %>%
group_by(capacity) %>%
top_n(1, loading) %>%
ungroup() %>%
select(-loading) %>%
mutate(factor_names = recode_factor(factor,
"MR1" = "Negative emotions",
"MR2" = "Cognition & control",
"MR3" = "Bodily sensations",
"MR4" = "Positive/social emotions"),
factor = factor(factor))
# make new dataframe
d_cat <- d_all %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(capacity, response, -c(subid, starts_with("target"))) %>%
left_join(factors_S2) %>%
left_join(d_demo %>%
select(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
mutate(subid = as.character(subid)))
Joining, by = "capacity"
Joining, by = "subid"
And get an average “score” for each of these factors for each participant:
d_cat_scored <- d_cat %>%
group_by(subid, parent,
target, target_num, target_ord,
factor, factor_names) %>%
summarise(score = mean(response, na.rm = T)) %>%
ungroup()
# bootstrapped means and 95% CIs
d_cat_scored_boot <- d_cat_scored %>%
group_by(target, target_num, target_ord, factor, factor_names) %>%
multi_boot_standard("score") %>%
ungroup()
ggplot(d_cat_scored,
aes(x = target_num, y = score, color = factor_names)) +
facet_grid(~ factor_names) +
geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_scored_boot, color = "black",
aes(y = mean, group = factor_names)) +
geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set2", guide = "none") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
scale_x_continuous(breaks = seq(0, 60, 12)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, untransformed",
x = "target age (months)")
ggplot(d_cat_scored,
aes(x = target_num, y = score, color = factor_names)) +
facet_grid(~ factor_names) +
geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_scored_boot, color = "black",
aes(y = mean, group = factor_names)) +
geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set2", guide = "none") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, square-root-transformed",
x = "age after square-root transformation (months)")
ggplot(d_cat_scored,
aes(x = target_ord, y = score, color = factor_names)) +
facet_grid(~ factor_names) +
geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_scored_boot, color = "black",
aes(y = mean, group = factor_names)) +
geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set2", guide = "none") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Age as an ordinal variable",
x = "age (ordinal)")
# contrasts(d_cat_scored$target_ord) <- contr.poly(13)
# contrasts(d_cat_scored$factor) <- contr.sum(4)
contrasts(d_cat$target_ord) <- contr.poly(13)
contrasts(d_cat$factor) <- contr.sum(4)
r1 <- lmer(response ~ factor * poly(target_num, 3) +
(1 + factor + target_num | subid),
d_cat %>%
mutate(target_num = scale(target_num, center = F)))
summary(r1)
Linear mixed model fit by REML ['lmerMod']
Formula: response ~ factor * poly(target_num, 3) + (1 + factor + target_num |
subid)
Data: d_cat %>% mutate(target_num = scale(target_num, center = F))
REML criterion at convergence: 536724.6
Scaled residuals:
Min 1Q Median 3Q Max
-5.5360 -0.4298 0.0225 0.5167 5.5499
Random effects:
Groups Name Variance Std.Dev. Corr
subid (Intercept) 200.49 14.159
factor1 156.13 12.495 0.65
factor2 145.83 12.076 -0.44 -0.59
factor3 69.75 8.352 -0.55 -0.34 -0.24
target_num 51.83 7.199 -0.64 -0.42 0.50 0.31
Residual 334.01 18.276
Number of obs: 61620, groups: subid, 237
Fixed effects:
Estimate Std. Error t value
(Intercept) 68.3739 0.7634 89.569
factor1 6.9512 0.8216 8.461
factor2 -38.7255 0.7947 -48.728
factor3 25.1428 0.5573 45.116
poly(target_num, 3)1 3241.2664 88.4095 36.662
poly(target_num, 3)2 -1237.5812 18.2759 -67.717
poly(target_num, 3)3 532.0777 18.2759 29.114
factor1:poly(target_num, 3)1 -544.0136 31.6548 -17.186
factor2:poly(target_num, 3)1 2670.8981 31.6548 84.376
factor3:poly(target_num, 3)1 -2550.8843 31.6548 -80.584
factor1:poly(target_num, 3)2 32.5828 31.6548 1.029
factor2:poly(target_num, 3)2 38.2604 31.6548 1.209
factor3:poly(target_num, 3)2 881.6516 31.6548 27.852
factor1:poly(target_num, 3)3 48.0128 31.6548 1.517
factor2:poly(target_num, 3)3 -391.9779 31.6548 -12.383
factor3:poly(target_num, 3)3 -354.1489 31.6548 -11.188
Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
This is weirdly hard to model… lots of alternative models didn’t converge. Also, those extremely high t-values make me a little suspicious, as do the identical standard errors for the last 9 terms. So, let’s take this with a bit of a grain of salt for now.
d_parent <- d_all %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(capacity, response, -c(subid, starts_with("target"))) %>%
left_join(d_demo %>%
select(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
mutate(subid = as.character(subid)))
Joining, by = "subid"
parent_counts <- d_demo %>%
distinct(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
count(parent) %>%
mutate(proportion = n/sum(n))
knitr::kable(parent_counts, digits = 3)
| parent | n | proportion |
|---|---|---|
| No | 146 | 0.616 |
| Yes | 88 | 0.371 |
| NA | 3 | 0.013 |
d_demo %>%
count(Parent, ChildrenOldestAge_collapse)
Given these numbers, I think our best bet it just to look at parents vs. non-parents, and not try to separate out parents of young children (too few!).
scores_S2 <- efa_S2$scores %>%
data.frame() %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(factor, score, -c(subid, starts_with("target"))) %>%
left_join(d_parent %>% distinct(subid, parent)) %>%
left_join(d_cat %>% distinct(factor, factor_names)) %>%
mutate(factor = factor(factor))
Joining, by = "subid"
Joining, by = "factor"
Column `factor` joining character vector and factor, coercing into character vector
# bootstrapped means and 95% CIs
d_parent_boot <- scores_S2 %>%
filter(!is.na(parent)) %>%
group_by(parent, target, target_num, target_ord, factor_names) %>%
multi_boot_standard("score") %>%
ungroup() %>%
data.frame()
scores_S2 %>%
filter(!is.na(parent)) %>%
ggplot(aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_parent_boot, aes(y = mean, group = parent)) +
geom_pointrange(data = d_parent_boot, fatten = 1.5, # note: too close to dodge
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = parent), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by parent status",
subtitle = "Exact age, untransformed",
color = "parent status: ",
x = "target age (months)", y = "score (NOTE: scales differ)")
scores_S2 %>%
filter(!is.na(parent)) %>%
ggplot(aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_parent_boot, aes(y = mean, group = parent),
position = position_dodge(width = 0.3)) +
geom_pointrange(data = d_parent_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
position = position_dodge(width = 0.3)) +
# geom_smooth(aes(group = parent), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by parent status",
subtitle = "Exact age, square-root transformed",
color = "parent status: ",
x = "age after square-root transformation (months)",
y = "score (NOTE: scales differ)")
scores_S2 %>%
filter(!is.na(parent)) %>%
ggplot(aes(x = target_ord, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_parent_boot, aes(y = mean, group = parent),
position = position_dodge(width = 0.5)) +
geom_pointrange(data = d_parent_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
position = position_dodge(width = 0.5)) +
# geom_smooth(aes(group = parent), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by parent status",
subtitle = "Age as ordinal variable",
color = "parent status: ",
x = "age (ordinal)",
y = "score (NOTE: scales differ)")
# bootstrapped means and 95% CIs
d_cat_parent_scored_boot <- d_cat_scored %>%
filter(!is.na(parent)) %>%
group_by(parent, target, target_num, target_ord, factor, factor_names) %>%
multi_boot_standard("score") %>%
ungroup()
ggplot(d_cat_parent_scored_boot,
aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_parent_scored_boot,
aes(y = mean, group = parent)) +
geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, untransformed",
x = "target age (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_parent_scored_boot,
aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_parent_scored_boot,
aes(y = mean, group = parent)) +
geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, square-root transformation",
x = "age after square-root transformation (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_parent_scored_boot,
aes(x = target_ord, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_parent_scored_boot,
aes(y = mean, group = parent)) +
geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Age as an ordinal variable",
x = "age (ordinal)",
y = "score (NOTE: scales differ)")
All of these visualizations suggest that, relative to non-parents (n = NA), parents (n = NA) tended to perceive greater abilities, in all domains except for bodily sensations. For both negative emotions and positive/social emotions, this difference is greatest in the mid-range of the target ages (12-24 months); for cognition & control, it’s greater at the older end of the target ages (e.g., 3-5 years). It doesn’t look like a huge effect, but it’s intruiging.
There are many ways that we could choose to model this - I’ll try out the factor scores first, adapting our primary regression models (not in this notebook).
contrasts(scores_S2$parent) <- contr.treatment(2, base = 1) # baseline: NON-parents
contrasts(scores_S2$factor) <- contr.sum(4)
r2 <- lmer(score ~ target_num * factor * parent
+ (target_num + factor | subid),
scores_S2 %>%
mutate(target_num = target_num/12))
summary(r2)
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ target_num * factor * parent + (target_num + factor |
subid)
Data: scores_S2 %>% mutate(target_num = target_num/12)
REML criterion at convergence: 22480.5
Scaled residuals:
Min 1Q Median 3Q Max
-9.3553 -0.4388 0.0524 0.5114 5.2974
Random effects:
Groups Name Variance Std.Dev. Corr
subid (Intercept) 0.3491 0.5908
target_num 0.0165 0.1284 -0.73
factor1 0.1947 0.4413 0.49 -0.40
factor2 0.2303 0.4799 -0.66 0.56 -0.40
factor3 0.2605 0.5104 0.24 0.00 -0.22 -0.61
Residual 0.2919 0.5403
Number of obs: 12168, groups: subid, 234
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.424399 0.049602 -8.556
target_num 0.279097 0.011334 24.624
factor1 0.078297 0.039260 1.994
factor2 -0.273518 0.042253 -6.473
factor3 0.253010 0.044634 5.669
parent2 0.070601 0.080884 0.873
target_num:factor1 -0.058596 0.006811 -8.603
target_num:factor2 0.190717 0.006811 28.001
target_num:factor3 -0.170985 0.006811 -25.104
target_num:parent2 0.003333 0.018482 0.180
factor1:parent2 0.005619 0.064020 0.088
factor2:parent2 -0.051089 0.068900 -0.741
factor3:parent2 0.031692 0.072783 0.435
target_num:factor1:parent2 -0.002939 0.011107 -0.265
target_num:factor2:parent2 0.053753 0.011107 4.840
target_num:factor3:parent2 -0.046204 0.011107 -4.160
Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
The last three coefficients here are where the action is: According to this model, the difference between parents and non-parents in perceived developmental trajectories (target_num:parent2, which is not significant collapsing across factors) is exaggerated in the domain of cognition and control (target_num:factor2:parent2) and diminished in the domain of bodily sensations (target_num:factor3:parent2).
Let’s try adding polynomial effects:
r3 <- lmer(score ~ poly(target_num, 3) * factor * parent
+ (poly(target_num, 1) + factor | subid),
scores_S2 %>%
mutate(target_num = target_num/12))
summary(r3)
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ poly(target_num, 3) * factor * parent + (poly(target_num,
1) + factor | subid)
Data: scores_S2 %>% mutate(target_num = target_num/12)
REML criterion at convergence: 19808.5
Scaled residuals:
Min 1Q Median 3Q Max
-10.3042 -0.4652 0.0017 0.5305 5.8317
Random effects:
Groups Name Variance Std.Dev. Corr
subid (Intercept) 0.2278 0.4773
poly(target_num, 1) 519.8012 22.7991 -0.51
factor1 0.1982 0.4452 0.45 -0.39
factor2 0.2338 0.4835 -0.60 0.54 -0.40
factor3 0.2640 0.5138 0.30 0.00 -0.22 -0.61
Residual 0.2316 0.4812
Number of obs: 12168, groups: subid, 234
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.030561 0.039884 -0.766
poly(target_num, 3)1 48.855288 1.983984 24.625
poly(target_num, 3)2 -19.550418 0.613116 -31.887
poly(target_num, 3)3 8.448995 0.613116 13.780
factor1 -0.004388 0.038065 -0.115
factor2 -0.004395 0.041145 -0.107
factor3 0.011731 0.043587 0.269
parent2 0.075304 0.065038 1.158
poly(target_num, 3)1:factor1 -10.257028 1.061949 -9.659
poly(target_num, 3)2:factor1 3.413024 1.061949 3.214
poly(target_num, 3)3:factor1 -1.785482 1.061949 -1.681
poly(target_num, 3)1:factor2 33.384462 1.061949 31.437
poly(target_num, 3)2:factor2 3.384910 1.061949 3.187
poly(target_num, 3)3:factor2 -6.433547 1.061949 -6.058
poly(target_num, 3)1:factor3 -29.930494 1.061949 -28.184
poly(target_num, 3)2:factor3 8.594647 1.061949 8.093
poly(target_num, 3)3:factor3 -3.251241 1.061949 -3.062
poly(target_num, 3)1:parent2 0.583455 3.235228 0.180
poly(target_num, 3)2:parent2 -2.031708 0.999792 -2.032
poly(target_num, 3)3:parent2 1.883539 0.999792 1.884
factor1:parent2 0.001472 0.062072 0.024
factor2:parent2 0.024763 0.067094 0.369
factor3:parent2 -0.033507 0.071076 -0.471
poly(target_num, 3)1:factor1:parent2 -0.514473 1.731690 -0.297
poly(target_num, 3)2:factor1:parent2 -0.041808 1.731690 -0.024
poly(target_num, 3)3:factor1:parent2 1.099314 1.731690 0.635
poly(target_num, 3)1:factor2:parent2 9.409337 1.731690 5.434
poly(target_num, 3)2:factor2:parent2 -2.361622 1.731690 -1.364
poly(target_num, 3)3:factor2:parent2 -1.057673 1.731690 -0.611
poly(target_num, 3)1:factor3:parent2 -8.087893 1.731690 -4.671
poly(target_num, 3)2:factor3:parent2 6.669873 1.731690 3.852
poly(target_num, 3)3:factor3:parent2 -4.830063 1.731690 -2.789
Correlation matrix not shown by default, as p = 32 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
A lot to sort through here, but here are some observations:
poly(target_num, 3)2:parent2 and poly(target_num, 3)3:parent2)poly(target_num, 3)1:factor2:parent2)poly(target_num, 3)1:factor3:parent2). I’m not sure how to interpret the differences in non-linearity here.You could imagine re-running any of these with a square-root transformation on target age, and/or with the ‘summary scores’ by category, but I’m not going to do that right now.
parent_counts_S1 <- demo_S1 %>%
distinct(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
count(parent) %>%
mutate(proportion = n/sum(n))
knitr::kable(parent_counts_S1, digits = 3)
| parent | n | proportion |
|---|---|---|
| No | 147 | 0.653 |
| Yes | 77 | 0.342 |
| NA | 1 | 0.004 |
scores_S1 <- efa_S1$scores %>%
data.frame() %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(factor, score, -c(subid, starts_with("target"))) %>%
left_join(demo_S1 %>% distinct(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
mutate(subid = as.character(subid))) %>%
mutate(factor_names = recode_factor(factor,
"MR4" = "Negative emotions (S1 F4)",
"MR1" = "Cognition & control (S1 F1)",
"MR2" = "Bodily sensations (S1 F2)",
"MR3" = "Positive/social emotions (S1 F3)"),
factor = factor(factor))
Joining, by = "subid"
# bootstrapped means and 95% CIs
d_parent_boot_S1 <- scores_S1 %>%
filter(!is.na(parent)) %>%
group_by(parent, target, target_num, target_ord, factor_names) %>%
multi_boot_standard("score") %>%
ungroup() %>%
data.frame()
scores_S1 %>%
filter(!is.na(parent)) %>%
ggplot(aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_parent_boot_S1,
position = position_dodge(width = 5),
aes(y = mean, group = parent)) +
geom_pointrange(data = d_parent_boot_S1,
position = position_dodge(width = 5),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = parent), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by parent status (STUDY 1)",
subtitle = "Exact age, untransformed",
color = "parent status: ",
x = "target age (months)", y = "score (NOTE: scales differ)")
I won’t bother to re-plot with different treatments of age, since there are only three target ages here.
# make new dataframe
d_cat_S1 <- d_all_S1 %>%
rename(subid_target = X) %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(capacity, response, -c(subid, starts_with("target"))) %>%
left_join(factors_S2) %>%
left_join(demo_S1 %>%
select(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
mutate(subid = as.character(subid))) %>%
filter(!is.na(factor))
Joining, by = "capacity"
Joining, by = "subid"
And get an average “score” for each of these factors for each participant:
d_cat_scored_S1 <- d_cat_S1 %>%
group_by(subid, parent,
target, target_num, target_ord,
factor, factor_names) %>%
summarise(score = mean(response, na.rm = T)) %>%
ungroup()
# bootstrapped means and 95% CIs
d_cat_parent_scored_boot_S1 <- d_cat_scored_S1 %>%
filter(!is.na(parent)) %>%
group_by(parent, target, target_num, target_ord, factor, factor_names) %>%
multi_boot_standard("score") %>%
ungroup()
ggplot(d_cat_parent_scored_boot_S1,
aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_parent_scored_boot_S1,
position = position_dodge(width = 5),
aes(y = mean, group = parent)) +
geom_pointrange(data = d_cat_parent_scored_boot_S1,
position = position_dodge(width = 5),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, untransformed",
x = "target age (months)",
y = "score (NOTE: scales differ)")
Looking quickly at these two visualizations of Study 1, I’d say that there’s a hint of the same patterns in the domains of negative emotions and positive social emotions (parents attributed more than non-parents, especially at 9 months), and in cognition/control (parents attributed more than non-parents at 5 years). So, generally consistent?
I won’t run regression analyses just now.